Relaxation and self-sustained oscillations in the time 
elapsed neuron network model 

Khashayar Pakdaman* Benoit Perthame^ Delphine Salort*"'- 

(N 

p ^ September 16, 2011 

<D 

in 

Abstract 

p ^ The izme elapsed model describes the firing activity of an homogenous assembly of neurons 

thanks to the distribution of times elapsed since the last discharge. It gives a mathematical de- 
. scription of the probability density of neurons structured by this time. In an earlier work, based on 

i^H generalized relative entropy methods, it is proved that for highly or weakly connected networks the 

model exhibits relaxation to the steady state and for moderately connected networks it is obtained 
numerical evidence of appearance of self-sustained periodic solutions. 
I— I Here, we go further and, using the particular form of the model, we quantify the regime where 

relaxation to a stationary state occurs in terms of the network connectivity. To introduce our 
^ methodology, we first consider the case where the neurons are not connected and we give a new 

statement showing that total asynchronous firing of neurons appears asymptotically. In a second 
T—i step, we consider the case with connections and give a low connectivity condition that still leads 

to asynchronous firing. Our low connectivity condition is somehow sharp because we can give an 
example, when this condition is not fulfilled, where synchronous rhythmic activity occurs. Indeed, 
On we are able to build several explicit families of periodic solutions. Our construction is fully nonlinear 

and the resynchronization of the neural activity in the network does not follow from bifurcation 
] analysis. It relies on an algebraically nonlinear boundary condition that occurs in the model. 

. . These analytic results are compared with numerical simulations under broader hypotheses and 

^ shown to be robust. 
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1 Introduction 



In nervous systems, neuronal circuits carry out tasks of information transmission and processing. 
During such functions, neuronal activity is evoked by incoming signals. However, in many assemblies, 
neuronal electrical activity is present and persists even in the absence of external signals. This form 
of activity that is not elicited by inputs is referred to as spontaneous activity (SA). That SA plays 
a pivotal role in functions such as respiratory rhythmogenesis that are and need to be maintained 
throughout the lifetime of organisms has been long recognized (see for instance [2j). Furthermore, for 
over the past twenty years, experimental evidence has been gathering for the widespread occurrence 
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of SA in different areas of nervous systems, such as the spinal cord, the retina, the cortex and the 
hippocampus, speciahy during various phases of development whereby SA mediates neuronal and 
network maturation (for reviews see |3][16]). In many instances, SA takes the form of recurring 
synchronous activity that emanates from a specific region and subpopulation of neurons and propagates 
throughout neuronal populations interacting through excitatory connections. In this work, our focus 
is on the mechanisms through which synchronous neuronal activity is initiated in tightly connected 
excitatory neuronal assemblies. 

Building on our previous mathematical analysis of dynamics of neuronal assemblies [17|, our main 
purpose is to use a mathematically tractable model of large neuronal populations incorporating a min- 
imal set of neuronal properties and construct explicitly periodic solutions of this model corresponding 
to synchronous oscillations reminiscent of SA. In the process, we also improve previous conditions 
ensuring asynchrony of assemblies. Our mathematical analysis and explicit construction are then 
compared with numerical simulations under broader hypotheses and shown to be robust. Our ap- 
proach highlights basic mechanisms lying at the core of synchronous rhythm generation in assemblies 
of interacting noisy excitable systems. 

That coupled oscillators can produce collective rhythms is well documented. Indeed, in such sys- 
tems, isolated units display periodic oscillations. Interactions between such units can entrain them 
into phase locking yielding global synchronous rhythms at the population level. The situation for 
coupled excitable systems is different. In this case, isolated units stabilize at a resting state and do 
not have intrinsic periodicity. It is in the presence of noise that they display occasional irregular 
discharges. Interactions between units can transform this irregular spontaneous activity of isolated 
excitable systems into a global coherent and regular rhythm. What makes this form of spontaneous 
regular synchronous activity of particular interest is that it emerges in a population of excitable rather 
than oscillating units. In other words, unlike synchrony in coupled oscillators, here synchrony and 
regularity of the global population are not direct consequences of intrinsic periodicity of the units. 

The emergence of such spontaneous synchronous oscillations in populations of noisy excitable neu- 
ronal models has been widely observed in the literature. One of the seldom cited earliest reports of 
this phenomenon is MacGregor and Palasek's study of randomly connected populations of neuromimes 
|13j . Other examples include the description of similar phenomena in models such as the active ro- 
tator [23! [Z], elementary neuron models [23], the FitzHugh-Nagumo [25], the Morris-Lecar [9], the 
Hindmarsch-Rose [5], the Hodgkin-Huxley |26j and even models implementing detailed biophysical 
properties |12j . In some of these references, the units are diffusively coupled, in others they interact 
through excitatory pulses; in some connectivity is all-to-all, while others deal with random networks. 
Our enumeration, which does not intend to be exhaustive, illustrates the ease with which assemblies of 
excitable units generate noise-induced synchronous regular activity, irrespective of model and network 
details. 

It is the widespread occurrence of the phenomenon in diverse network models that motivated us to 
design an elementary model with minimal assumptions shared by all excitable systems to gain better 
understanding of the way such rhythmic patterns arise. Previous explorations of the phenomenon are 
mainly numerical. Here we present the first example where there is an explicit construction of solutions 
for this problem. While our model was originally adjusted to experimental data |22j . in its present 
version, it contains only a bare skeleton of physiological neuronal properties. This reduction serves to 
emphasize and extract the key mechanisms at work. It shows that even with minimal assumptions it 
is possible to generate synchronous activity. Another advantage of the explicit construction is that it 
opens the door to study other forms of patterns that may arise in assemblies of interacting systems. 
We highlight this aspect by constructing other forms of periodic synchronous oscillations that had not 
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been reported previously. Using numerical simulations, we verify that the results are robust to the 
details of the model. 

From a mathematical standpoint, the time elapsed neuron population model belongs to the class 
of structured models that appear in a broad range of applications (see the reviews [131 EI])- The 
version considered in this paper possesses a threshold-like firing rate. This is similar to threshold like 
maturation and death rates in models of cell maturation, fish populations etc p!U]. In the same way 
as for these other models, this leads to a natural relation between our model and equations with state 
dependent delays. We take advantage of this property in the analysis of the dynamics of populations 
of unconnected neurons. However, for strongly interacting neuronal populations, no such relations 
exist and the analysis of the dynamics cannot be performed using this technique. In this respect, the 
periodic oscillations that are obtained in this work are also novel to structured populations and should 
be of broader interest in the field of structured models. 

This article is organized as follows. In section [2| we present the time elapsed neuron population 
model and discuss the underlying assumptions with respect to other related models. In the following 
section [Sj we consider the case of unconnected neuronal populations, henceforth referred to as neuronal 
ensemble. We prove that such a population tends to a state of total asynchronous neuronal activity. 
In our analysis , we use the underlying delay equation. This section also introduces tools that are used 
later in the paper. In section [4j we extend the conditions for asynchronous activity to interconnected 
populations of neurons satisfying a criterion on postdischarge recovery rate and the maximal refractory 
period. After these results, we focus on situations where the network can sustain synchrony. More 
precisely, in section [5} we build an explicit example of a network that possesses periodic solutions 
corresponding to rhythmic and synchronous activity in the network. Finally, in section [6], numerical 
results and variants of the model show the robustness of the results. We conclude the paper in section 
[7] by discussing our results, possible extensions and limitations. 

2 The time elapsed neuron population model 

The population model analyzed in this work encompasses a minimal set of properties reproducing the 
following aspects of neuronal behavior. Many neurons generate trains of stereotyped electrical pulses 
— referred to as spikes or action potentials — in response to incoming stimulations. Following each 
discharge, the neuron undergoes a period of refractoriness during which it is less responsive to inputs, 
before recovering its excitability 

One key aspect of neuronal coding is that the shape of the spike varies little. The main carrier 
of information is the discharge times or some statistics of the discharge times. Motivated by this, 
many neuronal and network models neglect the mechanisms underlying spike generation and instead 
represent neuronal dynamics in terms of discharge times. A general model of neuronal population 
dynamics based upon occurrence times of events was introduced in the pioneering work of Perkel 
\18\ \T9\ [20] . Other related variants are the spike response model or the integrate and fire models [6] . 
In this work, we consider a specific neuronal model that is inspired from our previous studies that 
aimed to reproduce SA in experimental preparations |22j . 

The fundamental assumption in our model is to describe the postdischarge recovery of the neuronal 
membranes through an instantaneous firing rate that depends on the time elapsed since the last 
discharge and the inputs by neurons. Such models have been used by others as well [6j. By adjusting 
this firing rate function, such models can reproduce experimental recordings of neuronal activity (see 
for instance [1]). Here, we adjust it for a different purpose, namely to construct a version in which 
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the occurrence of periodic oscillations is tractable. 

We consider a population of neurons described by the probability density n(s, t) of finding a neuron 
in 'state' s at time t where s represents the time elapsed since the last discharge. We assume here that 
the network is homogenous, that is to say that each neuron in the network has the same dynamic. 
Moreover, we assume that all the neurons are excitatory. According to |22| [23] a simple model of the 
dynamics is given by the age-structured equation 

\ N{t) := n{s = 0,t) = f^°^ p{s, X{t)) n{s, t)ds, 
completed with an initial probability density n^{s) that satisfies 

< n°(s) < 1, / n°(s)ds = 1. (2) 
Jo 

The coefficient p{s, X) represents the firing rate of neurons in the 'state s' and in an environment X 
resulting from the global neural activity. It is usually small (or even vanishes) during a rest phase that 
depends on X, and increases suddenly afterwards. The density of neurons undergoing a discharge at 
time t is denoted by N(t) and the boundary condition at s = means that the neuron re-enters the 
cycle at 'age s = 0' after firing. The interactions between the neurons are taken into account through 
the global neural activity at time t, X{t). In the theoretical part of this work we use 

X{t) = Nit) (3) 

which corresponds to instantaneous transmission of neuronal discharge. In the numerical part we also 
take into account synaptic integration, which is modeled by distributed delay function d, e.g.. 



X{t) = [ d{s)N{t- s)ds. 
Jo 

Here we also assume that the neurons can discharge earlier when the neural activity is high. We assume 
that the strength of connectivity is directly related to the sensitivity of the threshold of discharge with 
respect to the total neural activity. 

These modeling assumptions can be written on the (nonnegative) coefficients as 
d 

—p{s,x) > 0, pis,x) = for s G (0, s*(x)) p{s,x) > for s > s*{x), p{s,x) R{x), (4) 

d 

—p{s,x)>0, p{s,x)^^^l, < R{x) ^i^z^ 1. (5) 

Equations Q and ([s]) encompass the key biological properties incorporated into the model. The first 
is that following a discharge, a neuron looses its excitability, that is its propensity to discharge in 
response to a stimulation, and recovers it progressively in time. Furthermore, the mean firing rate of 
the neuron is a function of the stimulus level x, bounded by R{x). The second is that the probability of 
neuronal firing increases with the intensity of excitatory inputs, until some saturation level is reached. 
Here again, no further assumptions are made regarding the shape of p. 



4 



We will make constant use of two remarkable properties satisfied by solutions to ([T]) (see a 
conservation law which expresses the interpretation of n{t,-) as a probability density and an a priori 
bounds that reflects the normalization p{s, x) < 1, 

roo poo 

/ nis,t)ds= / n°(s)ds = l Vt > 0, (6) 
Jo Jo 

< n{s, t)<l, < N{t) < 1, \ft> 0, Vs > 0. (7) 

The theorems and the proofs of the paper are written for a specific choice of the function p which 
is piecewise constant. More precisely, p is given by 

pis,x) = I{s>a{x)}, for some a G Lip{R+), a'{x) < 0, 

(8) 

(7(0) =0-+, 0<a- < a{x) < a+, Vx > 0. 

Hence, at each time t, the neuron population undergoes two different dynamics as they are or not in 
their refractory state. 

• During the refractory period (s < a{N{t)), the population density satisfies the free transport 
equation. At each time t, the flux of neurons which enter in their refractory period is given by 
N{t) = n{s = 0,t). 

• After the refractory period, {s > a{N{t)), the solution satisfies the transport equation with 
an exponential decay due to the neurons discharge. At each time t, neurons enter this second 
stage where they begin to discharge. As we will see it later, the corresponding flux is given by 
n{a{N{t)),m-N'{t)a'{Nm. 

Therefore, at each time t, the function cr{N) defines the threshold time elapsed when the neurons begin 
to discharge. The derivative \a'\ measures the sensitivity of neurons to the global neural activity. Here 
our aim is to understand the link between \ a'\ and the appearance of rhythmic and synchronous activity 
in neural assemblies. Hence the quantity that we will essentially study in this article is not the function 
n{s,t), but the density of neurons which discharge at each time t, i.e., N{t). More precisely, we are 
going to focus on two dynamics. 

• The case where converges asymptotically to a constant: this dynamic corresponds to a total 
desynchronization of neuronal discharge in the network. 

• The case where is a periodic function: this dynamic corresponds to the appearance of rhythmic 
and synchronous activity, where the higher synchronization of neurons is given for the times t 
where A^ takes its bigger value. 

3 Dynamics of ensembles of disconnected neurons (^cr(x) = 0) 

To explain our method, we begin with the linear case, when neurons are firing independently from 
one another. We introduce a tool to prove desynchronization of the population. We first revisit and 
improve the results obtained in [T7] with using a mere norm and thus we remove the exponential 
weights on the initial data. Secondly, we propose a different analysis, by studying directly the firing 
rate function N{t). We prove that, as t tends to infinity, N(t) converges to a constant A^*, oscillating 
around this constant. 
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We assume that a is a small enough constant and thus the equation on n is linear given by 

f ^ + ^+p(«)n(.,t) = 0, 

\ N{t) := n{s = 0,t) = /o^'^p(s) n{s,t)ds, 

P{S) =I{s>a}, 0<a<l. (10) 

3.1 Long time behavior of the solution n{s,t) 

Theorem 1 With assumptions ^ and (10), the solution of the equation ^ converges with an ex- 
ponential rate to the stationary state 

A(s) = N*e--^oP^''^'^'\ with N* = ^— . (11) 

More precisely, there exist /i > and C > such that the following estimate holds 

r>+oo r+oo 

\n{s, t) - A{s)\ds < Ce-''* / |n(s, 0) - A{s)\ds. 
Jo 

Several variants of this result are known on bounded intervals (see [13] for instance). The interest 
here is to work on the half-line and in the functional space which is natural in the present context. 

Proof of Theorem [T]. 

Splitting of the solution in two terms. The idea of the proof consists in splitting the solution in 
two parts and set 

n° = n? + n(s, t) = ni{s, t) + n2{s, t). 

• The first initial term takes the value A^* on an interval strictly bigger than [0, a] and its support 
is not necessarily finite. 

• The second term has zero average and with bounded support, independent of (we can choose 
[0,3] for instance) and thus we can apply our earlier theory in [17J. 



To build this decomposition, we notice that the stationary solution A in (11) satisfies 

r+oo 

A{s) = N* on [0, a] and / A{s)ds = N* . 

J a 

Then, we split the initial data as follows 

with ^ 

5(nO)= r \\s)ds. 
Jo 

We will make use of properties of and 77-2, namely 

• /q^°° ni{s)ds = 1 and thus f^°° n2{s)ds = 0. 

• The support of 71-2 is contained in [0, 3] because a < 1 and a{nP) < 1. 
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• is equal to N* on an interval strictly bigger than [0,o"]. 

• n{s, t) - A{s) = n2{s, t) + [ni(s, t) - A{s)]. 

Time decay of n2- We recall the following Proposition proved in [I^, using an entropy method 

nam] 

Proposition 1 There exist fi > and a continuous function ip, uniformly hounded from below by a 
positive constant, such that for all bounded initial data with J^°° n^{s)ds = 0, the estimate holds: 

+00 r+oo 

|n(s,t)|V'(s)ds < e"^* / |n°(s)|V'(s)ds. 

^0 

We can apply Proposition [T] to the solution n2 and get 



+ 00 



P2 



is,t)\i;{s)ds <e-''^ / ?/;(s)|n^(s)|ds. 



Using the fact that the support on is contained in [0, 3], we obtain, with C 



SUP[0,3] ^(^) 

infa ip{s) 



+00 r+oo 

\n2is,t)\ds <Ce-^'^ \n%s)\ds. (12) 

Jo 



Estimate of ni — A. We begin with a Lemma that will be used on ni — A, 
Lemma 1 Let be such that /q^°° n^{s)ds = and = on [0, a]. Then, we have 

N{t) =n{s,t) = fors€[0,a], t>0. 
Proof of Lemma [l| Using the characteristics, we deduce that for all time t and s G [0, a] 

n{s, t) = N{t - s) where we extend for t G [-a, 0] by N{t) = 0. 
With this extension, the mass conservation gives 



N{t) + / N{t- s)ds = 0, Vt > 0. (13) 
Because N(t) = for t G [— u, 0], we conclude that for all r > 

sup |iV(OI < sup f \N{t- s)\ds <a sup \N{t)\ 

0<t<T 0<i<r70 0<t<T 



which implies (because a < 1) that = 0, which ends the proof of Lemma [T] □ 

To obtain the exponential decay in time of \\ni — A(t)\\j^i, we are reduced to proving the following 
Lemma (which we will apply to ni — j4 also) 



Lemma 2 Let nP G with f^°° n^{s)ds = be a function such that n^{s) = for s G [0, a] and let 
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n be the solution of the equation 



dn{s,t) dn{s,t) / f\ — n 

dt ^ ds ^ fH6,i; - u, ^^^^ 
N{t) := n{s = 0,t) = /+°° n{s,t)ds. 



Then the following estimate holds 



\n{s,t)\ds < e~* / (15) 
Jo 



Proof of Lemma [2| Since the function |n| is solution of 

dt\n\ + ds\n\ + \n\ = 0, 

using Lemma [T| we know that n{s,t) = for s G [0,ct]. Therefore we may integrate with respect to 
the variable s G [0, +oo[, we deduce that 



d 
dt 



+ 00 

\n\{s,t)ds = — I \n\{s,t)ds = — I \n\{s,t)ds 



JO 



from which we conclude the inequality (15) and Lemma [2] is proved. □ 
With Lemma [2] we conclude that 

Z' + OO /'OO 

/ \ni{s,t) - A{s)\ds <e-^ I \nl{s) - A{s)\ds. 
Jo Jo 

Combined with the estimate (12), the proof of Theorem [l] is complete. □ 
3.2 Oscillatory relaxation of the global activity 

As a preparation to the nonlinear case, we now explain another method to prove asymptotical total 
desynchronization. This method relies on a direct study of the function A^(-). 

We mention that, as we will see it later, we are in the simple situation where at each time t > a, 
the flux of neurons which enter in the area of discharge is given by 

nia,t)[l-^a'{Nit))]=N{t-a). 

More precisely the following Proposition holds: 

Proposition 2 With the assumptions of Theorem^ the solution N converges exponentially with rate 
In a to the constant N* = Moreover, N{t) oscillates around of N* ; for any interval I of size a 

3to G / such that N{to) = N* . 
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Proof of Proposition [2| Using the characteristics for equation ^ and mass conservation, we obtain 



N{t) + [ N{s)ds = 1, Vt > a, 



which we write as 

ft 



(16) 



N{t) -N* = [N* - N{s)]ds. 

Jt-a 

Using that < A'"* <l,0<A^<lin the right hand side, we obtain that 

\N{t)-N*\<a, yt>a. 
We may iterate the above argument and deduce that the following estimate holds 

|A^(t) - iV*| < ct", Vt>ncj, 
which ends the proof of the exponential convergence of N to N* because cr < 1. 
It is also easy to prove that the solution oscillates around N* . Let / = (toj + be an interval of 



size a. The relation (16) implies that 

N{to + a) + j^N{s)ds = l, 

which contradicts the equality A^* + //_^ N*ds = 1 if N{t) > N* or N{t) < N* on /. This completes 
the proof of Proposition [2j □ 

4 Asynchronous activity for connected neuronal assemblies 7^ 0) 

We come back to the nonlinear situation when the threshold depends on the global neural activity to 
study relaxation properties. In |17| it was proved that both for weakly and highly connected networks, 
the model undergoes relaxation to the unique steady state. But the precise conditions in there are 
difficult to follow and are certainly not optimal. Our purpose here is to come to precise (and somehow 
sharp) conditions leading to global relaxation. We introduce a different strategy, following the method 
of section 13.21 

Our first condition is that for some < m < 1 we have 

< -ct'(x) < m < 1. (17) 
It implies that for any probability density < n(-) < 1, there is a unique N such that 

/■oo 

N = n{s)ds, 0<N <1, 

Ju{N) 

because i— /j^^y) ^{s)ds is a contraction and we can apply the Banach-Picard Theorem. In particular 
the boundary condition ([T]) defines a unique N{t), while the periodic solutions we build later correspond 
to two possible choices of N, hence discontinuities. 
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4.1 Some observations 

The conservation law ^ can be written 

1 = N{t) + [ n{s,t)ds. 
Jo 

The following lemma ensures that N is smooth enough for our purpose 



(18) 



Lemma 3 Assume that condition (17) holds. Then, the function N is Lipschitz continuous. 



Proof of Lemma [sj Consider two times ti, t2- The relation (18) gives for i = 1,2 



1 = N{ti) + [ 
Jo 



n{s, ti)ds. 



hence 



N{ti)-N{t2)= [n{s,t2) -n{s,ti)]ds+ n{s,t2)ds. 

Jo JcT{N{ti)) 



Using the fact that < n < 1 and condition (17), we find that 



n{s, t2)ds 

/<T(7V(tl)) 

On the other hand, using the PDE on n, we have 



<m\Niti)-N{t2)\. 



(19) 



/c 



<T(7V(tl)) 





(s, ta)- n{s, ti)]ds = f^^ dtn{s, t)dsdt 

= & r^^^*^^^ d,n{s, t) + Pis, N{t))n{s, t)dsdt 



112 JO 

= fo'^''^'"^^ P('^ N{t))n{s, t)dsdt + Hi [n{cj{N{t)),t) - N{t)]dt. 

We obtain that there exists a constant C such that 

r(T{N(ti)) 

10 



[n{s,t2) - n{s,ti)]ds 



< C\ti-t2\. 



(20) 



Combining estimates (19) and (20) and the fact that m < 1, we deduce that there exists a constant 
C{m) such that 

\N{ti)-N{t2)\<C{m)\ti-t2\ 
which ends the proof of Lemma [3j □ 



We may now differentiate a.e. the relation ( 18 ) and find 

dN{t) 



dt 



■ [1 + a'{N{t))n{aiN{t)), t)] = -Nit) + n(a(7V(t)), t) . 



(21) 



This equation is much more complicated to use than in the case when a is constant. It can be 
interpreted as follows. At each time t, 



10 



• the boundary condition tells us that the flux of neurons which enter in their refractory state is given 
by N{t). 

• the flux of neurons which enter in the state of possible discharge is given by 

n{a{Nmt) [l-^a'{N{t))]. 



This quantity is not as easy to control as in the linear case. For instance it can take nonpositive values 
if the velocity of the threshold given by '^^^'^ a' {N) gets bigger than the chronological speed. 

4.2 A condition for asymptotic total desynchronization 



In this section, using assumption (17), we give a sufficient condition on the function a and the maximal 
refractory duration 

= sup (7(2;) = cr(0) 

R+ 

that implies asymptotical total desynchronization of the neural activity. 



Theorem 2 Assume (17) and that, with N <1 the unique steady state defined in (26) below, 

a{x) < cr+ < 1 - mN, 
Then, we have '^^i^^ cr' {N (t)) < 1 and with an exponential rate. 



(22) 



lim N(t) = N. 

In sectionjsj we give explicit examples where the equality —^p-cr' {N{t)) = 1 holds and N{t) can be 
periodic. 

Proof of Theorem [2| Assumption (17) and the a priori estimate < n{s,t) < 1 imply that 

1 > g{t) := 1 + a'{N{t))n{a{N{t)),t) > 1 - m > 0. 



Hence, using the relation (21), we conclude 

dN{t) _ 
dt ~ 



-N{t)+n{a{N{t)),t) 

git) 



We decompose the end of the proof as follows 

In a first step, we prove that condition (17) implies that the threshold velocity is always strictly less 

(23) 



than 1 as soon t > , i.e., the following estimate holds 



^^^^^^-a'iNit)) <m, \ft> a+, 



dt 



a property (23) that allows us to simplify the equation on A^. 

• In a second step, this simplified equation on can be used to show that A^ converges, when t — )• +00. 
This is based on an iterative process similar to that used in section [S] 

As a first step we prove the 
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Lemma 4 With condition (u7[), the estimate (23) holds. 



Proof of Lemma [4j Using the relation (21) on A^, we deduce that, with the notation n{t) 
n{a{N{t)),t)) > 0, we have (still because a' < 0) 



dNjt) 
dt 



a 



'{N(t)) 



a'[N{t)){n(t)-N{t)\ 



1 + 



l+(T'(Af(t)) n(t) 

l-a'(jVft))jV(t) 



l+<T'(Af(t)) n{t) 

< 1 + (—1 + m) = m. 



Next, we simplify the equation on N using that the threshold velocity is strictly less than 1. 



□ 



Lemma 5 Assume that estimate {23) holds. Then, we have 

ft 



Nit) + / N{s)ds = 1 Vt > CJ+. 

Jt-a(N{t)) 

Proof of Lemma [5| Using the characteristics, we obtain that 

n{s,t) = N{t - s)e-^oP('''^(^+''-'^)'^'' for s < a{N{t)). 
Take s € [0, a{N{t))] and let us prove that 

n{s,t) = N{t - s). 

To do this, it enough to show that in the above integral p{s' , N{t + s' — s)) = 0, that is 

Vs' G [0, s) we have s' < a{N{s' + t- .s)). 

We set 

f{s') ■.= s'-a{N{s' + t-s)). 



We have /(s) < because s < a{N{t)). Moreover, from (23), we deduce that 

dN{s' + t-s) , 



f{s') = 1 



dt 



-a'{N{s' + t- s) > 0, 



which proves (25). We conclude the proof of Lemma [s] using (18) which now reads 



1 = N{t) + / n(s, t)ds = N{t) + / 



N{t - s)ds. 



(24) 



(25) 



□ 



We can now conclude the proof of Theorem [2] We notice there is a unique steady state solution 
< < 1, that is a solution to 

N{l + a{N)) = l. (26) 
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This is because the mapping N 1 — Na{N) is a contraction as soon asO<o"(-)<l and — 1 < 
f'CO < and thus it has a unique fixed point by the Banach-Picard theorem. 
Then, for all t > o""*" we rewrite (24) as 

N{t)-N + 
This proves that for t > 



t-a{N(t)) 



[N{s) - N]ds + N[a{N{t)) - a{N)] = 0. 



\N{t) -N\- N\a{N{t)) - (j{N)\ < 



t-a{N{t)) 



\N{s)-N\ds, 



\N{t)-N < 



1 - mN t^aiNit))<s<t 



sup |iV(s)-iV|. 



Since < N,N{-) < 1, we find \N{t) - N\ < jf^. 
We may iterate and for t > na^ 

\N{t)-N\ < 



a ' 



1 - mN 



Using condition (22), this ends the proof of Theorem 2 



□ 



5 Synchronous periodic oscillations in neuronal assemblies 

In our previous study [T7], we have presented numerical examples of periodic solutions corresponding 
to synchronous activity that were obtained in numerical simulations. Here, for a particular class of 
functions a{-) we construct analytically such periodic solutions. To do this, we impose a priori that, 
on a given open interval /, or a union of two intervals, neurons do not leave their refractory state, 
namely 

|a(iV(i))=l VtG/. (27) 

This can only last until saturation occurs, leading suddenly to a massive discharge of neurons charac- 
terized by a discontinuity of N{t). 

For a > 0, we define the two functions 

< N-{a) := < N+{a) := — - < 1, (28) 

V ^ 2e° - 1 ^ ^ 2e" - 1 ^ ' 



and we choose the Lipschitz continuous discharge threshold a as 
a{x) 



' 2a on [0,iV^(a)], 

2a-ln(x) + ln(iV-(a)) on [N-{a),N+{a)\, (29) 
a on [A^+(a),oo). 



Motivated by numerical results that we show in section [6j we are going to build three distinct classes 
of periodic solutions. The first two are simpler than the third, but they seem unstable because only 
the third class is accessible by our numerical simulations. As we will see, there are several common 
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features in our constructions. One of them is to postulate a form of N{t) and solve the linear transport 
equation 

f ^ + ^ + Pis, N{t)) nis,t)=0, tGR,s> 0, ^^^^ 
\ n{s = 0,t)= N{t). 

Because it is periodic in t, there is no need of initial data and the method of characteristics gives a 
solution. The three classes depend on suitable choices of N{t). 

5.1 Periodic solutions with one discontinuity 

The first class is characterized by a periodic activity N{t) with a single discontinuity per period. 

Theorem 3 For a > 0, define the discontinuous periodic function N{t) of period a as 

N{t) = N+{a)e-^ for tel := (0,a). 



Then the solution to the linear renewal equation (30) is also a solution of the equation ([ijj. In other 
words, we have for all t gM 

/■oo fOO 

N{t) = / n[s,t)ds and / n{s,t)ds = I. 

Ja{N(t)) Jo 

Some comments about this Theorem are in order 
• Since N{a~^) = N~{a), N(t) S (N~{a),N~^{a)) is a sawtooth function; it decreases exponentially 
from N~^(a) to (a). When the discharge rate of neurons is at its minimal value N~(a), a peak of 
discharge occurs in the network which is characterized by a jump of the function from N~{a) to 
N~^{a). Consequently the precise shape of a outside the interval [N~{a), N'^{a)] does not play a role 
here. 



We can compute that, for these solutions, the assumption (17) does not hold. Indeed 



sup W'{N{t))\ = l/N-{a) = 2e" - 1 > 1. 



• The bigger a is, the bigger is the difference N'^{a) — N (a). In particular, A^+(0) — N (0) = and 

lim A^+(a) = - and lim N~{a) = 0. 

o— ^+oo 2 a—>--\-co 

The synchronization phenomenon is more evident for large values of a. Hence we can hope that sus- 
tained oscillations still occur for more general models when a is large enough. For instance, numerical 
simulations indicate that this choice of the function a also yields periodic solutions when the global 
neural activity is defined with a delay as in |17j . 

Proof of Theorem [sj To begin with, we notice that (27) holds true with our choice. Then we argue 
in two steps. First we establish a form of the boundary condition, then the integral condition. 
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First step. We notice that, for t £ I, the discontinuities of n{s,t) lies on {s = ka + t}. Therefore, 



taking into account (p7|), we calculate, for t S / 

A 

dt 



= n{0,t) =N{t) = -N'{t). 
We may now write, by periodicity, that for all t € M, 



a{Nit)) 

n{s,t)ds + N{t) ■.= Qi = l. (31) 







Indeed, in order to evaluate the constant Qi, we can use the method of characteristics and compute 
this quantity at t = 0+ 

l-cr{N{0+)) i-cr{N+{a)) pa 

/ n{s,0)ds= / N{-s)ds = N+{a) / e"-"ds = Af+(a)(l - e^"). 

Jo Jo Jo 

After evaluating this expression, we arrive at 

r-<7(7V(0+)) 



n(s, 0)ds = 1 - iV+(a), Qi = 1. 
Second step. Finally, integrating equation (|30[), we compute for t G M (the second equality uses 







(3Tf) 

d 
dt 



poo nOG roo 

/ n{s,t)ds = N{t) — / n{s,t)ds = l— / n{s,t)ds. 

Jo Ja{N{t)) Jo 



For a periodic function this proves that n(s, t)ds = 1 and the proof of Theorem [s] is complete. □ 
5.2 Periodic solutions with two discontinuities 

It is possible to build a family of more elaborate periodic solutions by allowing for two discontinuities 
of N(t) during a period. The construction uses an additional parameter 

p G (0,a), Np{a) = iV_(a) G (iV_(a), 7V+(a)), (32) 

and the single root ^{p) £]p + a, 2a] to the equation (7(7) = where 

g(j) = e-"[e"-P - l][e^-" - 7 + 1] + e~"[-2a + (a + p)e"-P] + e^"" - 1. (33) 

Theorem 4 For p £ (0, a), define the ^{p) periodic function by 

7V+(a) e"* for t € [0,a], 

Ar(t) = ,( iVp(a) e-*+" /or t G (a,p + a], 

7V+(q) e-*+°[t-7(p) + eT(?')-"] /or t G (p + a, 7(p)). 



Then the solution to the linear renewal equation (30) is also a solution of the equation ([ijj. In other 
words, we have for all t £M 

fOO POO 

N{t) = / n{s,t)ds and / n{s,t)ds = l. 
Ja(N{t)) Jo 
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The single peak solution in section 5.1 corresponds to the limiting case p = 0, 7(p) = a. 



Proof. As in the proof of Theorem [3| we have for t E (0, a) and t G {a,p + a 

— / n{s,t)ds = N'a'{N{t))n{a{N{t)),t)+ dt 

Jo Jo 

j-a{N{t)) 

= n{a{N{t)),t) - dsn{s,t)ds 



= n{0,t) = N{t) = -N'{t). 
For t £ (a + p,^{p)), we have N{t) > N^{a). Therefore a{N{t)) = a and we calculate 

n{s,t)ds = N'a'{N{t))n{a{N{t)),t)+ / dtn{s,t)ds 



dt 



a{N{t)) 

dsn{s, t)ds 



= n{0,t) -n{a,t) = N{t) - N{t - a) 

= iV+(a)e-*+"[t - 7(p) + e'^P^] - iV+(a)e-*+" = -N'{t). 

By periodicity, we conclude that there is a constant m{p) such that for all t G M 

n{s,t)ds + N{t) = m{p). 



Furthermore, following the second step in the proof of Theorem [3j we have for all t > 

»+oo 

n{s, t)ds = m{p). 



It remains to check that m{p) = 1 and it is enough to prove that 

/ n(s,0+)ds + A^(0) = 1. 

Jo 

That is to say 

— - / N{s)ds + eP-'' = 1. 
-(a) Ja+p 



K 
We have 

nip) 



r Nis)ds = e-"[e"-P - l][e^(P)-" - -f{p) + 1] + e-"[-2a + (a +p)e"-P]. 

-(") Ja+p 



With the definition in (33), we deduce that we must have g{'^{p)) = 0. 

It remains to prove that there exists a unique value 7(1?) G (a + p,2a) such that g{'^{p)) = 0. To 
prove this, we observe that 

g{a+p) = e-"[e"-f + p - a - 1] < 0, 

and that 

g{2a) = e-"e"-P[e" +p-a + l]-2 + e?'"" - e"" > 0. 
Since g{-) is increasing, this completes the proof of Theorem |4| □ 
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5.3 Periodic solutions with two discontinuities and a flat state 



We are now ready to present the third class of periodic solutions which seems to be those observed 
numerically. We use again the notations (|32l) and we need the function 



f{6,y) := e°([-2a + l + e" + 5]e- 



-2a 



e"") + e-"[5-ln(y)-2a]+y-l. 



(34) 



Theorem 5 For p £ (0,a) and 6{p) G [p + a,2a], let Y{p) be such that f{6{p),Y{p)) = (we will see 
that 5{p) exists), let N he the 2a periodic function given by 



' N+{a) e- 



Nit) 



for t G [0,a], 

for t £ {a,p + a], 

7V_(a) for t G [p + a, 6{p)], 

N+la) e-*+"[t - 2a + e"] for t £ {5{p),2a). 



Np{a) e-*+" 



Then the solution to the linear renewal equation (30) is also a solution of the equation In other 
words, we have for all t gM 

poo roc 

N{t) = / n{s,t)ds and / n{s,t)ds = 1. 
Ja(N{t)) Jo 

The double peak solution in section 5.2 corresponds to 6 = p + a and Y{p) = e^"" (no flat state). 



Proof. As before, we can write for t G (0, a) and t & {a,p + a 
dt 



aiNit)) 



dsn{s, t)ds 



n{s,t)ds = N'a'{N{t))n{a{N{t)),t) - 
= n{0,t) = N{t) = -N'{t). 
For t G {a+p,6{p)), N{t) = iV_(a) and so cr{N{t)) = 2a on {a + p,6{p)). Hence for t G {a + p,6{p)), 



di 



n{s,t)ds = N'a'{N{t))n{a{N{t)),t) 

aiN{t)) 



aiN(t)) 



dsn{s, t)ds 



dsn{s, t)ds 

= n(0, t) - n{2a, t) = N{t) - N{t - 2a) = = -N'{t). 
For t G {6{p),2a), we have N{t) > N+{a). We deduce that cr{N{t)) = a and that 



dt 



n{s,t)ds = N'a'{N{t))n{a{N{t)),t) 

a{N{t)) 



dgU^s, t)ds 



dsn{s, t)ds 

= n{0,t) -n{a,t) = N{t) - N{t- a) 

= iV+(a)e-*+"[t -2a + e"] - 7V+(a)e-*+" = -N'{t). 

Again, we conclude by periodicity that there is a constant m{p) such that for all t G M 

ra(7V(t)) 
10 



f 

Jo 



n{s,t)ds + N{t) = m{p). 



(35) 
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To prove that for all t > 

»+oo 

n{s, t)ds = m{P), 



we use the equahty ( 35 ) , and we have 



p+oo r+oo p+oo 

/ n{s,t)ds = N{t) — / n{s,t)ds = m{p) — I n{s,t)ds. 

Jo Ja{N{t)) Jo 

By periodicity, we deduce that for all t > 

X) 

n(s, t)ds = m{p). 

Let us assume that the following lemma holds and let us prove that m{p) = 1. 



Lemma 6 Recalling the function f in (34), for allY{p) € [e 1], there exists 5{p) G [p + a,2a] such 
that, 

f{6{p),Y{p)) = 0. 
To prove that m{P) = 1, it is enough to check that 

r<T{N(0+)) 

n(s,0+)ds + Af(0) = 1, 

N{s)ds = 1 - N+{a). (36) 
N{s)ds = Np{a)P - N-{a) 



that is 

r-2a 



We have 

fP+a 



hence, by setting Y{p) = e^ we rewrite estimate (36) as 

r2a 



We have 

r-2Q 



\-- I N{s)ds + e-'^[6{p) - ln(y(p)) - 2a] + Y{p) = 1. 
-(a) Js{p) 

N[s)ds = iV+(a)e" f [-2a + 1 + e" + 5{p)]e-^^P^ - e'^" - e"") . 



'S(p) 

Hence the relation ( 36 ) is equivalent to the following equality 

e"f [-2a + 1 + + 6{p)]e-^^P'^ - e'^" - e"") + e-°[(5(p) - HY{p)) - 2a] + Y{p) = 1. 



As 6{p) is chosen such that f{S{p),Y{p)) = 0, this conclude the proof of Theorem [5] assuming that 
Lemma [6)iolds. □ 

Proof of Lemma [6[ We argue by continuity in 5 and consider the endpoints. On the one hand, we 
have 

/(a + p, Y{p)) = ^ [1 + e" + ln(y (p))] - e"" - e'^" + Y{p) - 1. 
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We are going to prove that 

/(a + p,y(p))>/i(e-")>0. 
with the function h : [e~", 1] — )• M given by 



To do so, we compute 



h{Y) := ^[1 + e° + ln(y)] - e"" - e'^" + Y-1. 



V(y) = i + .-"™O>0. 



We deduce the inequahty. 
On the other hand, 

/(2a, Y{p)) = -e-" ln(y(p)) + Y{p) - 1 < 0, 
which concludes the proof of Lemma [6j □ 

6 Simulations 

This section deals with the numerical solution to the nonlinear renewal equation ^ with the transi- 



tion from refractory state given by the piecewise smooth function (28) that has been used to derive 
analytical periodic solutions in section [sj The two questions we address here are 1) which of these 
periodic solutions can we observe numerically, 2) what is the effect of smoothing terms that describe 
memory effects such as synaptic integration. In the present, synaptic integration occurs through a 
delay between the actual firing rate N{t) defined as in ([T]) and its input in the firing rates coefficients 
X{t) in section [l| Rather than the simple rule (|3]), one introduces synaptic integration either through 
a convolution, as it was done in [17], or through a differential equation which is easier for numerical 
simulations 

X^+X{t) = N{t). (37) 

Here A > is interpreted as the relaxation time associated with the synaptic integration through a 
passive membrane. Our purpose is to see the effect of this parameter in the stability of the periodic 
solutions built in section [H 



6.1 Numerical methods 

We compute the numerical solution with a standard upwind scheme. The theory of such finite volume 
schemes can be found in [H [8] and in the references therein. Computationally it gives a very simple 
algorithm where the solution n{s, t) is approximated by ~ ^ Ss''^i^2 ^^^^^ with grid parameters 
defined by As = At (this means that the Courant-Friedrichs-Lewy (CFL in short) number equal to 
one, as it is usual for age-structured equations. We also use the notation = (i + 1/2) As and 

= kAt. The nonlinearities specific to the model ([l|, ^ are treated as follows. The firing rate term 
p(s, X(t))n(s, t) is treated implicitly (and thus with unconditional stability) which allows us to keep 
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Figure 1: (Instantaneous transmission) Total neural activity N{t) as computed with the numerical 
scheme (38)~([40)) with the data in section [sj Left: a = .2. Right: a = .5. The continuous lines give 
the values and N^. This is a numerical solution compatible with our construction of solutions in 



section 5.3 with p > 0. 






16 18 20 22 24 26 28 30 32 34 36 



Figure 2: (Instantaneous transmission) Total neural activity N{t) (left) and solution n(s) (right) 
computed as in Figure [T] with a = 3. This numerical simulation is compatible with the exact solutions 



in section 5.3 with p > 0. 
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Figure 3: (Instantaneous transmission) As in Figure [2] with the initial data n^{s) = e ^. Comparison 
shows that various periodic solutions are attractive depending on the initial data. Here we obtain 



periodic solutions as in section 5.3 with p = 



the maximal time step needed for accuracy in the transport term. Then one finds the discrete values 
(n^)i<j</^^ iteratively on the time variable thanks to 



n 



k+l 



n: 



k ._ 



l<i< I, 



M, 



(38) 



1=1 



with p^=I 



{iAs>a{X'')} 



. It has to be understood that the label i = stands for the boundary value 



(s = at the continuous level) and thus the formula (38). The semi-implicit scheme does not preserve 
the fundamental conservation law (|6]) indicating we work with probabilities. Therefore one introduces 
a second order correction step (which in practice seems to have very little effect) and set 



n. 



fc+i 



When there is no delay we just use the idendity 



X 



k+l 



When time delay is included according to (37), we use an explicit Euler scheme and write 

At, 



X 



k+l 



x^ii 



(39) 



(40) 



(41) 



Numerical solutions are computed with 1000 grid points per unit length (we have tested that this 
allows to reach numerical convergence). 



6.2 Numerical periodic solutions 

The numerical results are displayed in Figures [T] and [2] for times such that the periodic regime is 
established. It seems the period is precisely 2a. The numerical solution is therefore not the one peak 
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Figure 4: (Time delay) Total neural activity N{t) (left) and solution n(s) (right) computed as in 
Figure [T] with a = 3. But time delay is included with A = .1. 



periodic solution built in section 5.1 with period a. It is not either the periodic solution of period less 



than 2a build in section 5.2 This fact led us to build the third class in section [5^ which precisely fits 



with the figures. The flat state can be observed even in figure [2] (zooming if necessary). 



The effect of synaptic integration is to stabilize the dynamics and periodic solutions do not appear 
when A is too large. When periodic solutions are still computed, the effect is to smooth out the 
solution n{s) away from the discontinuity which is still present. 



6.3 Transition from desynchronization to sustained activity 

We have also included a parameter J that describes the connectivity of the network. This can be 
implemented by simply using the threshold a{JN{t)) in place of a{N(t)) and 

p{s,N) = I{s>cr{JN))}- 

The numerical scheme shows a rapid transition at a critical value J*; for J < J* there is desyn- 
chronization and for J > J* discontinuous periodic solutions of large amplitudes appear similar to 
those depicted previously. This confirms that the periodic solutions do not result from a supercritical 
Hopf-bifurcation. 



7 Conclusion 

The 'time elapsed' model for neural networks is a model that uses the probability density of neurons 
structured by the time elapsed since the last discharge and that has been introduced in [IT] and is 
based on the stochastic simulations in [22l[23]. Following the study in [17\, both desynchronization or 
periodic solution may occur depending on the nonlinearity. This nonlinearity takes into account the 
connectivity of the network through a modulation of the refractory period depending upon the total 
activity of the network. 

Here we have quantified the nonlinearities that lead to total desynchronization. We arrive to a 
condition which is somehow sharp because we can build periodic solutions when the previous condition 
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is not fulfilled. Numerical simulations show that the periodic solutions are numerous and depend on 
the initial data and the analytical solutions built in section [5] are not always the stable ones observed 
numerically. 

We conjecture that the periodic solutions are always discontinuous (both the total network activity 
N{-) and the probability distribution of neurons n) in accordance with the numerical computations 
shown in section |6] and with the theoretical construction in section [5j 



References 

[1] Berry M. J., Meister M. Refractoriness and neural precision. J. Neurosci. 18(1998) 2200- 2211. 

[2] Bianchi, A. L., Denavit-Saubi M. and Champagnat J. Central control of breathing in mammals: 
neuronal circuitry, membrane properties, and neurotransmitters. Physiol. Rev. 75 (1995)1-45. 

[3] Blankenship A. G. and Feller, M. B. Mechanisms underlying spontaneous patterned activity in 
developing neural circuits. Nature Reviews Neuroscience 11 (2010) 18-29. 

[4] Bouchut F. Nonlinear stability of finite volume methods for hyperbolic conservation laws and well- 
balanced schemes for sources. Series Frontiers in Mathematics, Birkhauser Verlag, Basel (2004). 

[5] Du-Qu W., Xiao-Shu L. Coherence resonance and noise-induced synchronization in Hindmarsh 
Rose neural network with different topologies. Commun. Theor. Phys. 48 (2007) 759-762. 

[6] Gerstner W. and Kistler W., Spiking Neuron Models (Single Neurons, Populations, Plasticity. 
Cambridge Univ. Press, Cambridge (2002). 

[7] Giacomin G., Pakdaman K., Pellegrin X. and Poquet C. Transitions and active rotator systems: 
invariant hyperbolic manifold approach. Submitted. 

[8] Godlewski E. and Raviart P.-A. Numerical approximation of hyperbolic systems of conservation 
laws. Applied Math. Sciences 118, Springer- Verlag New- York (1996). 

[9] Han S. K., Yim T. G., Postnov D. E. and Sosnovtseva O. V. Interacting coherence resonance 
oscillators. Phys. Rev. Lett. 83 (1999) 1771-1774. 

[10] Hbid M. L., Snchez E. and Bravo de la Parra R. State-dependent delays associated to threshold 
phenomena in structured population dynamics. Math. Models Methods Appl. Sci. 17 (2007), no. 
6, 877-900. 

[11] Kandel E. R., Schwartz J. H. and Jessell T. M. Principles of Neural Science, 4th ed. (2000) 
McGraw-Hill, New York. 

[12] Kosmidis E.K. , Pierrefiche O., Vibert J.-F. Respiratory-like rhythmic activity can be produced 
by an excitatory network of non-pacemaker neuron models. J. Neurophysiol. 92 (2004) 686-699. 

[13] MacGregor R. J. and Palasek R. L. Computer simulation of rhythmic oscillations in neuron pools. 
Biological Cybernetics 16 (1972) 79-86. 

[14] Metz J. A. J. and Diekmann O. The dynamics of physiologically structured populations. LN in 
biomathematics 68 Springer- Verlag (1986). 



23 



[15] Michel P., Mischler S. and Perthame B. General relative entropy inequality: an illustration on 
growth models. J. Math. Pures Appl. 84(9):1235-1260 (2005). 

[16] Moody W. J. and Bosma M. M. Ion channel development, spontaneous activity, and activity 
dependent development in nerve and muscle cells. Physiol. Rev. 85 (2005) 883-941. 

[17] Pakdaman K., Perthame B. and Salort D. Dynamics of a structured neuron population. Nonlin- 
earity 23 (2010) 55-75. 

[18] Perkel D. H. A computer program for simulating a network of interacting neurons. I. Organization 
and physiological assumptions. Comput Biomed Res. 9:31-43 (1976). 

[19] Perkel D. H. A computer program for simulating a network of interacting neurons. III. applica- 
tions. Comput Biomed Res. 9:67-74 (1976). 

[20] Pcrkcl D. H. and Smith M. S. A computer program for simulating a network of interacting neurons. 
II. Programming aspects. Comput Biomed Res. 9:45-66 (1976). 

[21] Perthame B., Transport equations in biology. Series 'Frontiers in Mathematics', Birkhauser (2007). 

[22] Pham J., Pakdaman K., Champagnat J. and Vibert J.-F. Activity in sparsely connected excitatory 
neural networks: effect of connectivity. Neural Networks 11 (1998) 415-434. 

[23] Pham J., Pakdaman K. and Vibert J.-F. Noise-induced coherent oscillations in randomly con- 
nected neural networks. Phys. Rev. E 58(3) (1998) 415-434. 

[24] Shinomoto S. and Kuramoto Y. Phase transitions in active rotator systems. Prog. Theor. Phys. 
75 (1986) 1105-1110. 

[25] Toral R., Mirasso C, Gunton J. D. System size coherence resonance in coupled FitzHugh-Nagumo 
models. Europhysics Letters 61 (2003) 162-167. 

[26] Wang Y., Chik D. T. W. and Wang Z. D. Coherence resonance and noise-induced synchronization 
in globally coupled Hodgkin-Huxley neurons Phys. Rev. E 61 (2000) 740-746. 



24 



